/*
 * exportVTK2.cpp
 *
 *  Created on: 12/07/2011
 *      Author: rogsoares
 */
#include "exportVTK.h"


void printScalarField(ofstream &, pMesh, MapSCDBL::const_iterator);
void printVTKHeader(ofstream&, pMesh);
void printVerticesCoordenates(ofstream&, pMesh);
void printElementConnectivities(ofstream&, pMesh);
void printElementConnectivities(ofstream&, pEntity, const int &);
void printCellTypeList(ofstream&, pMesh );
void printPressure(ofstream&, pMesh);
void printScalarField(ofstream&, pMesh, string);

void exportSolutionToVTK(pMesh theMesh, string filename, const MapSCDBL &mapData){
	
  cout << __LINE__ <<":\t" << __FILE__ << endl;
  
	// open file
	ofstream fid;

	cout << __LINE__ <<":\t" << __FILE__ << endl;
	
	fid.open(filename.c_str());
      
	cout << __LINE__ <<":\t" << __FILE__ << endl;
	
	// print data to file
	fid << "# vtk DataFile Version 2.0\n";
	
	cout << __LINE__ <<":\t" << __FILE__ << endl;
	
	fid << "Two phases flow simulation\n";
	fid << "ASCII\nDATASET UNSTRUCTURED_GRID\n";
	fid << "POINTS " << M_numVertices(theMesh) << " float\n";

	cout << __LINE__ <<":\t" << __FILE__ << endl;
	
	printVerticesCoordenates(fid,theMesh);
	
	cout << __LINE__ <<":\t" << __FILE__ << endl;
	
	printElementConnectivities(fid,theMesh);
	printCellTypeList(fid,theMesh);

	// loop over all properties strings to be printed in the VTK file
	fid << "\nPOINT_DATA "<< M_numVertices(theMesh) << endl;
	MapSCDBL::const_iterator iter;
	for(iter = mapData.begin(); iter != mapData.end(); iter++) printScalarField(fid,theMesh,iter);
	fid.close();
}

void printScalarField(ofstream &fid, pMesh theMesh, MapSCDBL::const_iterator iter){
	fid << "SCALARS " << iter->first << " float 1\n";
	fid << "LOOKUP_TABLE default\n";
	pEntity node;
	VIter vit = M_vertexIter(theMesh);
	while( (node = VIter_next(vit)) ){
		int i = EN_id(node);
		fid << iter->second[i] << endl;
	}
	fid<<"\n";
	VIter_delete(vit);
}

void exportSolutionToVTK(pMesh theMesh, string filename, const std::list<string> &propStrLst){
	// open file
	ofstream fid;
	fid.open(filename.c_str());

	// print data to file
	fid << "# vtk DataFile Version 2.0\n";
	fid << "Two phases flow simulation\n";
	fid << "ASCII\nDATASET UNSTRUCTURED_GRID\n";
	fid << "POINTS " << M_numVertices(theMesh) << " float\n";

	printVerticesCoordenates(fid,theMesh);
	printElementConnectivities(fid,theMesh);
	printCellTypeList(fid,theMesh);

	// loop over all properties strings to be printed in the VTK file
	fid << "\nPOINT_DATA "<< M_numVertices(theMesh) << endl;
	for(std::list<string>::const_iterator iter = propStrLst.begin(); iter != propStrLst.end(); iter++ )
		printScalarField(fid,theMesh,*iter);

	fid.close();
}

void printScalarField(ofstream &fid, pMesh theMesh, string field){
	fid << "SCALARS " << field << " float 1\n";
	fid << "LOOKUP_TABLE default\n";
	pEntity node;
	double val;
	VIter vit = M_vertexIter(theMesh);
	while( (node = VIter_next(vit)) ){
		EN_getDataDbl(node,MD_lookupMeshDataId(field.c_str()),&val);
		fid << val << endl;
	}
	fid<<"\n";
	VIter_delete(vit);
}


// print vertex coordenates and transfer to each one all computed values
// to be printed later
// = = = = = = = = = = = = = = = =  = = = = = = = = = = = = = = = = = =
void printVerticesCoordenates(ofstream &fid, pMesh theMesh){
	pEntity e;
	int count = 0;
	VIter vit = M_vertexIter(theMesh);
	while( (e = VIter_next(vit)) ) {
		double coord[3] = {.0, .0, .0};
		V_coord(e,coord);
		for (int i=0; i<3; i++) fid << coord[i] << " ";
		fid << endl;
		EN_attachDataInt(e,MD_lookupMeshDataId("mLN"),count++);
	}
	VIter_delete(vit);
}

// print elements connectivities
// = = = = = = = = = = = = = = = =  = = = = = = = = = = = = = = = = = =
void printElementConnectivities(ofstream &fid, pMesh theMesh){
	pEntity e;
	const int dim = theMesh->getDim();
	const int numElements = (dim==2)?M_numFaces(theMesh):M_numRegions(theMesh);
	fid << "\nCELLS " << numElements << " " << (dim+2)*numElements << endl;

	if (dim==2)
	{
		FIter fit = M_faceIter(theMesh);
		while( (e = FIter_next(fit)) )
			printElementConnectivities(fid,e,dim);
		FIter_delete(fit);
	}
	else
	{
		RIter rit = M_regionIter(theMesh);
		while( (e = RIter_next(rit)))
			printElementConnectivities(fid,e,dim);
		RIter_delete(rit);
	}
}

void printElementConnectivities(ofstream &fid, pEntity e, const int &dim){
	int mappedLNodes[4];
	std::vector<pVertex> vertex;
	M_GetVertices(e,vertex);

	fid << dim+1 << " ";
	for(int i=0; i<dim+1; i++){
		EN_getDataInt(vertex[i],MD_lookupMeshDataId("mLN"),&mappedLNodes[i]);
		fid << mappedLNodes[i] << " ";
	}
	fid << endl;
	vertex.clear();
}

void printCellTypeList(ofstream &fid, pMesh theMesh){
	const int dim = theMesh->getDim();
	const int numElements = (dim==2)?M_numFaces(theMesh):M_numRegions(theMesh);
	const int celltype = (dim==2)?5:10;
	fid << "\nCELL_TYPES " << numElements << endl;
	for(int i=0; i<numElements; i++) fid << celltype << endl;
}

void printPressure(ofstream &fid, pMesh theMesh){
	// print scalar field
	// = = = = = = = = = = = = = = = =  = = = = = = = = = = = = = = = = = =
	fid << "\nPOINT_DATA "<< M_numVertices(theMesh) << endl;
	fid << "SCALARS scalar_field float 1\n";
	fid << "LOOKUP_TABLE default\n";
	pEntity node;
	double val;
	VIter vit = M_vertexIter(theMesh);
	while( (node = VIter_next(vit)) ){
		EN_getDataDbl(node, MD_lookupMeshDataId("interpScalar"), &val);
		//val = pPPData->getPressure(node);
		fid << val << endl;
	}
	VIter_delete(vit);
}
